Australian  Government 
Department  of  Defence 

Defence  Science  and 
Technology  Organisation 


Quadratic  Optimisation  with  One  Quadratic  Equality 

Constrai  nt 


H  atem  H  mam 

Electronic  Warfare  and  Radar  Division 

Defence  Science  and  Technology  Organisation 

DSTO-TR-2416 


ABSTRACT  (U) 

This  report  presents  a  theoretical  framework  for  minimising  a  quadratic  objective  function 
subject  to  a  quadratic  equality  constraint.  The  first  part  of  the  report  gives  a  detailed  algorithm 
which  computes  the  global  minimiser  without  calling  special  nonlinear  optimisation  solvers. 
The  second  part  of  the  report  shows  how  the  developed  theory  can  be  applied  to  solve  the 
time  of  arrival  geolocation  problem. 


RELEASE  LIMITATION 


A  pproved  for  public  release 


Published  by 

Electronic  Warfare  and  Radar  Division 

DSTO  Defence  Science  and  Technology  Organisation 

PO  Box  1500 

Edinburgh  South  Australia  5111  A  ustralia 

Telephone:  (08)7389  5555 
Fax:  (08)  7389  6567 

©  Commonwealth  of  A  ustralia  2010 

AR-014-772 

June2010 


APPROVED  FOR  PUBLIC  RELEASE 


Quadratic  Optimisation  with  One  Quadratic 
Equality  Constrai  nt 


Executi  ve  Su  mmary 

The  theoretical  work  presented  in  this  report  is  motivated  by  the  need  to  solve  many 
defence  application  problems  such  as  the  time  of  arrival  geolocation  problem.  The 
mathematical  tools  needed  to  solve  such  a  localisation  problem,  aredeveloped  in  detail 
in  the  first  part  of  the  report  (Part  I)  and  are  based  on  quadratic  optimization  with  one 
quadratic  equality  constraint.  The  main  contribution  of  this  work  is  the  development  of 
an  algorithm  which  provides  a  step-by-step  procedure  to  solve  the  problem  and 
address  solution  feasibility  and  uniqueness  issues.  This  algorithm  relies  heavily  on 
linear  algebraic  transformations  and  optimality  condition  properties  to  efficiently  and 
exactly  determine  the  problem  minimiser.  The  localisation  of  an  emitter  source  or  a 
receiver  based  on  time  of  arrival  (TOA)  measurements  is  a  demonstration  example 
given  in  the  second  part  of  the  report  (Part  II),  which  illustrates  how  the  developed 
theory  is  used. 
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1.  Introduction 


The  current  work  was  motivated  by  the  need  to  solve  nonconvex  optimisation  problems 
encountered  in  many  engineering  applications.  This  report  focuses  on  developing  an 
efficient  algorithm  to  exactly  solve  the  problem  of  minimising  a  multivariate  quadratic 
cost  function  subject  to  a  single  quadratic  equality  constraint.  The  class  of  problems  we 
would  like  to  solve  is 


minimize  z  P  z  +  2p  z  +  n 

Z 

subject  to  zT Q z+2qT z  +  v  =  0 


(1.1) 


where  the  variable  z  e  R"  and  Q  e  S" ,  (Q  ^  0) .  The  parameters  p  and  q  are  vectors  of  R" 
while  u  and  v  are  real  numbers.  The  objective  Hessian  matrix  P  is  assumed  positive 
definite.  An  application  where  optimisation  of  this  type  arises  is  global  positioning  of  a 
receiver  using  pseudo-range  measurements.  Section  6  of  the  report  briefly  demonstrates 
how  this  and  related  applications  fit  the  general  framework  of  (1.1)  by  applying  the 
algorithm  developed  in  this  report  to  a  time-of-arrival  (TOA)  localisation  problem. 

Variants  of  problem  (1.1)  appeared  in  the  mathematical  literature  [3]  [4]  but  not  with  the 
intent  to  give  a  detailed  algorithm  to  solve  it.  In  fact,  as  we  shall  see  later,  advanced  linear 
algebra  techniques  are  used  to  transform  (1.1)  into  a  form  where  the  problem  can  be  solved 
efficiently  and  exactly.  Fundamental  results  on  the  first  and  second  order  optimality 
conditions  [3] [4]  help  locate  the  global  minimiser  of  (1.1). 

This  report  is  arranged  as  follows.  Section  2  reformulates  (1.1)  using  advanced  linear 
algebra  techniques.  The  purpose  is  to  render  optimisation  problem  (1.1)  more  tractable. 
Section  3  presents  the  notations  used  in  this  work.  It  also  defines  key  matrices,  eigenvalues 
and  vectors  that  shape  the  behaviour  of  the  optimisation  problem.  Section  4  addresses 
problem  feasibility  issues  related  to  the  constraint  function.  Section  5  presents  optimality 
conditions  and  provides  an  algorithm  that  solves  (1.1)  efficiently  and  exactly.  Section  6 
forms  the  second  part  of  the  report  and  gives  an  illustrative  geolocation  example  showing 
how  the  developed  constrained  optimisation  algorithm  is  applied.  Section  7  provides 
concluding  remarks. 


2.  Problem  Reformulation 


To  solve  problem  (1.1),  we  transform  (1.1)  into  a  more  tractable  problem  as  follows. 
Because  P  is  positive  definite,  then  there  exists  an  invertible  matrix,  S,  which 
simultaneously  satisfies  SPST  =In  and  SQST  =QD  (diagonal)  [7],  The  matrix  I n  is  the 
nth  dimensional  identity  matrix.  Using  this  matrix  decomposition,  problem  (1.1)  can  be 
equivalently  re-expressed  as 
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minimize  -v  A,  -v  +  ^  -v  + w 

y 

subject  to  _y+2f r  v  +  v  =  0 


(1.2) 


where  the  parameters  of  (1.1)  are  related  to  those  of  (1.2)  through  z  =  STy ,  b  =  Sp  and 
f  =  Sq .  Using  the  substitution,  x  =  y  +  b,  this  equation  can  in  turn,  be  transformed  into 


T  T 

minimize  A  «  x  +  r 

x 

subject  to  xtQd  x+2 ctx  +  5  =  0 


(1.2') 


where  r  =  u  —  bTb,  c  =  f  -QDb  and  s  =  v  +  bTQDb-2fTb  .  Given  that  the  minimiser  of 
(1.2')  is  not  affected  by  the  value  of  r,  we  may  select  r  =  0 ,  and  therefore  the  basic 
optimisation  problem  to  be  solved  reduces  to 


minimize  x  x 

X 

subject  to  h(x )  =  xTQD  x+2 cTx  +  5  =  0 


(1.3) 


The  remainder  of  the  report  focuses  on  finding  globally  optimal  solutions  for  problem  (1.1) 
through  solving  (1.3).  In  particular  if  x *  is  a  global  minimum  solution  of  (1.3),  then  the 
minimiser  of  (1.1)  is  deduced  as  z*  =  ST  (x*  -  Sp) .  Note  that  (1.3)  can  also  be  viewed  as  the 
problem  of  finding  the  closest  point  on  a  n-dimensional  conic  to  the  origin.  Furthermore  it 
should  be  noted  that  if  c  =  0  in  (1.3),  then  only  square  terms  exist  in  both  the  constraint 
and  objective  functions  of  (1.3).  Solving  such  a  problem  becomes  straightforward  as  (1.3) 
can  be  transformed  into  a  linear  program  (LP).  In  the  remaining  part  of  the  report,  we 
consider  the  non-trivial  problem  where  c  +  0 . 


3.  Some  Useful  Terminologies  and  Notations 

In  this  section  we  introduce  a  number  of  notations,  which  prove  useful  in  solving  (1.3)  or 
(1.1).  These  notations  are  related  to  QD  and  C.  First  note  that  if  QD  is  negative  semidefinite, 
then  by  multiplying  the  constraint  of  (1.3)  by  -1,  the  constraint  remains  unchanged  but  QD 
becomes  positive  semidefinite.  Flence  in  the  rest  of  the  paper  we  won't  consider  negative 
semidefinite  QD  .  By  convention  the  diagonal  matrix,  QD  ,  is  represented  as 
diag(//1;- ••,//,,  jUM  0---0)  where  the  m  nonzero  entries,  //,  ,  are  arranged  in  a 

descending  order  as  >  /u2  >•••>//,>  0  >  juM  >■■•>  jum  .  Note  that  because  QD  is 
congruent  to  Q  (i.e.  SQST  =  QD  with  S  invertible)  then  by  Sylvester's  law  of  inertia  [7], 
Qd  has  the  same  number  of  positive,  negative  and  zero  eigenvalues  as  Q .  We  also  define 
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the  largest  positive  and  smallest  negative  eigenvalues  of  QD  as  //max  >  0  and  jumm  <  0 .  As 
an  example  if 


4 

"0“ 

4 

0 

-  1 

and  c  = 

4 

0 

0 

0 

2 

then  1  =  2,  m  =  3,  //niax  =  4  and  jumm  =  - 1  .  We  also  denote  by  cmax  and  cmin  the 
subvectors  of  C,  which  correspond  respectively  to  the  extreme  eigenvalues  //rnax  and 


u  of  On  .  Hence  c  = 

•  min  max 


and  cmm  =  [4] .  Finally  Qn  is  singular  in  this  example  and 


therefore  we  let  c,  be  the  component  vector  of  C,  associated  with  the  zero  eigenvalue  of 


Qd  .  This  leads  to  c.  = 


0 

2 


in  this  example. 


4.  Problem  Feasibility 


In  this  section  we  introduce  a  theorem  detailing  the  conditions  under  which  (1.3)  is 
feasible. 

Theorem  1. 


1. 

2. 

3. 

4. 


If  Qd  is  indefinite,  then  the  feasibility  set  of  (1.3)  is  non-empty 

If  Qd  >  0 ,  singular  and  c,  ^  0,  then  the  feasibility  set  of  (1.3)  is  non-empty 


If  Qd  >  0 ,  then  the  feasibility  set  of  (1.3)  is  non-empty  if  and  only  if  s  <  V  — 

.  i  Mi 

If  Qd  >  0  and  c,  =  0,  then  the  feasibility  set  of  (1.3)  is  non-empty  if  and  only  if 

m  2 

s<Y  — 

M  Mi 


Proof. 


The  constraint  equation  in  (1.3)  can  be  rearranged  in  the  form 

/  m  n 

Kx)  =  Yj  Vixf  +  X  Mixf  +  2  X  cixi  +  s 

i= 1  /=/+!  i=m+ 1 
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where  /  corresponds  to  the  number  of  positive  eigenvalues  (counting  multiplicities).  The 
second  sum  of  squares  corresponds  to  the  negative  eigenvalues  of  Qn  .  The  third  term  (if 
any)  corresponds  to  the  zero  eigenvalues  of  Qn  . 

1.  If  Qd  is  indefinite  then  clearly  h(x)  =  0  has  always  real  solutions  and  hence  feasible 
points  always  exist. 

2.  If  Qd  >  0  and  singular,  then  there  are  no  negative  eigenvalues  (i.e.  m  =  / )  and  therefore 
no  second  quadratic  term  appears  in  h(x)  .  If  c:  ^  0  then  there  exists  at  least  one  non-zero 
entry,  ck ,  in  the  vector,  c: .  The  associated  variable,  xk ,  has  no  corresponding  square  term 
in  h(x) .  This  means  that  by  varying  xk  and  fixing  the  remaining  variables,  h(x)  can  be 
made  arbitrarily  large  (positive)  or  small  (negative)  and  therefore  by  virtue  of  continuity  of 
h(x)  ,  solutions  always  exist  to  h(x)  =  0  .  Feasible  points  always  exist  under  these 
conditions. 


3.  If  Qd  >  0 ,  then  the  minimum  of  h(x)  is  rj  =  s  -  V  — .  Hence  feasible  points  exist  if  and 

/=1  Mi 


only  if  r/  <  0 ,  or  in  other  words  s  <  V  —  . 

/=i  Mi 


4.  Similarly  if  QD>0  and  c_  =  0  ,  then  the  minimum  of  h(x)  is  r]  =  s  —  and 

/=i  Mi 


ZC  ■ 

— — . 

/= 1  Mi 


Theorem  1  gives  the  conditions  under  which  (1.3)  is  feasible.  These,  as  explained  below, 
also  include  uninteresting  cases  where  constraint  feasibility  is  restricted  to  a  single  point 

n  £.2 

(i.e.  when  Qn  >  0  and  s  =  /— )  or  to  a  subspace  of  reduced  dimension  (i.e.  when 

/=i  Mi 


Qd>0  singular,  c. 


The  minimisation  solution  for  these  two  special  cases  turns  out  to  be  trivial  and  unique  as 
shown  below. 


"  c2 

1st  case:  QD>  0  and  s  =  V  — 

/=!  Mi 


The  only  feasible  point  is  clearly,  x  =  —QDlc ,  which  is  the  global  minimiser  of  (1.3). 
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m  ^2 

2nd  case:  QD  >  0  singular,  s  =  V  —  and  c.  =  0 

M  A, 

First  we  define  the  matrix  2d«z  1°  be  the  submatrix  of  Qn  containing  all  nonzero  diagonal 
entries  and  then  we  denote  by  cnz  the  subvector  of  C  corresponding  to  QDnz .  If  we  partition 


x  as  x  =  nz  ,  then  clearly  the  point  xnz  =  -  QDircnz  satisfies  the  constraint  equation  with 

•v,  J 

x_  arbitrary.  It  is  straightforward  to  see  that  the  minimiser  of  (1.3)  then  becomes 

*  r  Q  Dnz^  nz  1 


To  rule  out  these  uninteresting  rare  cases,  we  also  require  that  [3]  [4] 

inf  h(x )  <  0  <  sup  h(x)  (1.4) 

R"  R" 

Condition  (1.4)  always  holds  for  scenarios  1  and  2  of  Theorem  1,  but  holds  in  scenarios  3 
and  4  if  and  only  if  the  inequalities  are  replaced  with  strict  inequalities.  In  the  rest  of  the 
report,  we  assume  that  problem  feasibility  holds  in  the  sense  of  (1.4). 


5.  Finding  the  Global  M  inimum  of  (1.1) 

Before  solving  (1.1)  or  equivalently  (1.3)  we  need  to  be  aware  that  the  solution  behaviour  of 
(1.3)  is  affected  by  c  v  and  c  .  .  Recall  that  these  are  the  subvectors  of  C  associated 
respectively  with  the  largest  positive  and  smallest  negative  eigenvalues  of  Q .  Appendix  A 
illustrate  how  cmax  =  0  or  cmin  =  0  leads  to  solving  a  reduced  optimisation  problem,  often 
with  no  unique  solution.  This  observation  has  motivated  the  introduction  of  the  next 
assumption. 

Assumption  1. 

Referring  to  problem  formulation  (1.3), 

if  Qn  >  0  ,  then  we  assume  that  cmax  ^  0  , 

if  Qd  is  indefinite,  then  we  assume  that  cmax  ^  0  and  cmin  ^  0  . 

This  assumption  gives  a  sufficient  condition  for  the  global  minimum  to  be  unique  (a  single 
point  in  R" )  and  helps  to  avoid  the  so-called  "hard  case"  problem  in  [6],  Going  back  to  the 
example  given  in  the  first  part  of  this  section,  QD  has  both  positive  and  negative 
eigenvalues.  We  clearly  have  cmin  ^  0  but  we  don't  have  cmax  ^  0  .  Hence  Assumption  1  is 
not  satisfied  in  this  example. 
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Next  we  define  the  Lagrangian  of  problem  (1.3)  as 
L(x,  /l)  =  xTIn  x  -  A(xTQD  x  +  2 cT x  +  s ) 

An  optimal  solution  must  satisfy  the  following  first-order  and  constraint  conditions 


(ln  ~  2 ,Qd  )x  -  Ac  =  0 
xtQd  x  +  2 cT x  +  5  =  0 


(1.5) 


From  [3]  [4]  [5],  the  global  minimum  must  also  satisfy  the  second  order  optimality 
condition 


I,-2Qd>0  (1.6) 

If  A  is  an  optimal  Lagrange  multiplier  satisfying  (1.5)  then  from  (1.6)  we  must  have 
Dx  =  In  —  A  Qd  >  0 .  The  next  proposition,  helps  us  narrow  down  the  region  where  the 

optimal  Lagrange  multiplier,  A  ,  resides. 


Proposition  1. 


U  QD>0  ,  then  A*  <  Amax  =  — . 

Amax 

1 

If  Qd  is  indefinite,  then - =  Amin  <A<  Amax 

Amin 


1 


A 


max 


Proof. 


The  proof  mainly  relies  on  the  second  order  optimality  condition  In  —  A  QD  >0  at  the 
global  minimum.  This  condition  leads  to  1  —  A  //,  >  0  for  all  nonzero  eigenvalues  of  QD  . 
Hence  the  Lagrange  multiplier  region  is  given  as  the  intersection  of  all  intervals  generated 
by  each  eigenvalue  inequality  ( 1  —  X /ut  >  0  ).  This  necessarily  means  that  the  Lagrange 
multiplier  region  is  an  interval  delimited,  from  one  side  or  both,  by  the  inverse  of  the 
largest  positive  or  smallest  negative  eigenvalues  of  QD  as  indicated  in  the  proposition.  □ 

This  proposition  shows  that  the  optimal  Lagrange  region  is  an  interval  of  finite  or  infinite 
length.  Let's  denote  this  interval  by  A  =  \Amin  A]mx  ]  ,  where  it  should  be  understood  that 
Arnin  can  be  finite  or  infinite  and  A  -  <  0  <  Am„  . 

null  min  niaA 


The  first-order  optimality  condition  in  (1.5)  reads  as 
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(I-XQd)x  =  Xc  (L7) 

Clearly  given  Assumption  1  and  Proposition  1,  it  necessarily  follows  from  (1.7)  that  the 
optimal  Lagrange  parameter,  X  ,  cannot  be  Xmj„  or  Xmax .  In  other  words  X  must  reside 
within  the  open  interval  A  =  (Xmm  Xmax )  and  D ..  is  invertible. 

From  (1.7),  the  globally  optimal  vector  x  can  be  obtained  as  a  function  of  X  as 
x  =  XD^c  (1.8) 

Plugging  (1.8)  into  the  constraint  xT QD  x+2 cTx  +  5  =  0  of  (1.5)  leads  to 


K(X)  =  X2ctD/QdD:1c  +  2XcTD-lc  +  5  =  0 


(1.9) 


which  reduces  to 


m 

K(X)  =  Y 


■  +  2- 


Xcf 


(i-M)2  (i-M) 


+  2/t  c2  +  5  =  0 


(1.10) 


i=m+ 1 


or  equivalently 


t t  jUi^-XjUi)  ti£  1 


2X^c2  +7  =  0 


(1.11) 


where  77  =  5  —  '5+  — . 

.=1  A, 

The  linear  term  in  (1.11)  is  present  if  and  only  if  (9fl  is  singular  and  c_  +  0 . 

The  next  proposition  examines  the  number  of  solutions  of  the  secular  equation  (1.11) 
within  the  open  interval  A . 


Proposition  2. 


If 

1.  Qd  is  indefinite,  or 

2.  Qd>0  singular  and  cz  +  0 ,  or 

3.  Qd  >  0  singular,  c_  =  0  and  77  <  0 ,  or 

4.  go  >  0  and  rj  <0 , 


then  K(X)  admits  a  unique  finite  solution,  X,  in  A .  Furthermore  the  constraint  function 
K(X )  is  increasing  with  X  . 
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Proof. 

First  we  prove  the  existence  of  a  solution  A  in  A  to  K (A)  =  0 .  We  will  carry  out  this 
analysis  by  distinguishing  whether  Qr:  is  positive,  indefinite,  singular  and  so  on. 


If  Qd  is  indefinite,  then  from  Assumption  1,  we  have  cmax  +  0  and  cmm  +  0  in  (1. 11).  Also 
following  Proposition  1,  we  have  Amm  =  — - —  and  =  — - —  .  Therefore 


lim  K(A ) 

=  -oo  and 

lim  K(A ) 

L^mm  J 

_A— >/lmax 

Amin 
=  +00 


An 


If  Qd  >  0  but  singular,  then  we  have  m  <  n  .  From  Assumption  1,  we  also  have  cmax  +  0  in 
(1.11).  Using  Proposition  1  it  necessarily  follows  that  Amm  =  — oo  and  Amdx  =  — - —  .  Hence 


~] 

f-oo  if  c,  +  0 

- 

lim  K(A)\ 

=  \  "  and 

lim  K(X) 

L^min  J 

[rj  <  0  if  c,  =  0 

_X— >/lmax 

If  Qd  >  0,  then  we  have  m  =  n  and  cmax  +  0  in  (1.11).  From  Proposition  1,  Amin  =  -oo  and 


/Lmax  = - )  and  it  follows  that 

Amax 


lim  K(A) 

=  r/  <  0  and 

lim  K(A ) 

L^^min  J 

Clearly  we  conclude  that  in  all  cases 


lim  K(A) 

X 

lim  K(A) 

J 

_A— >/^max 

<  0 ,  and  therefore  by 


continuity  of  K(A) ,  it  necessarily  follows  that  there  must  be  at  least  one  finite  solution  A , 
in  A,  satisfying  K(A)  =  0 . 


Furthermore  the  derivative  of  K(A)  with  respect  to  A,  is - =  2 which  is 

dA  ,=i  (1  -A/if) 

clearly  positive  for  all  A  in  A,  implying  that  K(A)  is  increasing  and  therefore  K(A) 
admits  only  one  finite  root,  A  ,  in  A . 

□ 

Proposition  2  shows  that  the  function  K(A)  is  well  behaved.  Figure  1.  shows  one  example 
of  how  K(A)  might  look  like  when  the  constraint  function,  QD,  is  indefinite.  It  is  clear  that 
irrespective  of  the  size  n  of  problem  (1.3)  or  (1.1),  efficient  zero  finding  algorithms  may  be 
applied  and  are  guaranteed  to  converge.  Note  also  that  A  =  0  is  always  within  A,  which 
may  be  used  as  a  starting  point  for  a  Newton-based  zero-finding  algorithm.  Furthermore  if 
Qd  is  indefinite  (as  in  Figure  1.)  then  a  bisection  algorithm  may  be  used  to  find  the  zero  of 
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K(A) .  The  starting  interval  for  the  bisection  algorithm  is  [Amm  +s  Amax  —  fjcz  A,  where 
£  >  0 ,  small  enough  to  be  able  to  evaluate  K(A)  at  the  first  iteration.  To  speed  up 
convergence,  a  switch  to  a  faster  zero  finding  algorithm  may  be  enabled  after  a  few 
iterations  of  the  bisection  algorithm.  If  however,  QD  is  positive  semidefinite  then  K(A)  is 


convex  within  A  =  (-00  Amax )  (because  the  second  derivative 


d~K 

dA2 


a  -Aii,y 


is 


positive  in  A)  and  the  zero-finding  algorithm  may  be  tailored  to  take  advantage  of  the 
convex  property. 


Figure  1:  An  illustrative  example  of  how  K(A )  might  look.  The  optimal  Lagrange  multiplier  is 
within  A  =  (-1/2  1/3) 
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The  algorithm  to  solve  (1.1)  is  summarised  below: 

1.  Ensure  that  Hessian  matrix,  P ,  of  the  objective  function  of  (1.1)  is  positive  definite. 
Otherwise  stop 

2.  Transform  problem  (1.1)  into  problem  (1.3)  using  standard  simultaneous  matrix 
decomposition.  If  QD  <  0 ,  multiply  the  constraint  equation  by  -1.  The  new  QD 
becomes  -  Q[: ,  which  is  positive  semidefinite. 

3.  Examine  the  feasibility  of  (1.3)  according  to  Theorem  1  and  property  (1.4).  If  (1.3)  is 

n  c2 

not  feasible  then  stop.  If  QD  >  0  and  s  =/\  — ,  or  if  QD>  0  singular,  c.  =  0  and 

,=i  A 

m  ^2 

s  =  — ,  then  carry  out  the  special  solving  procedure  described  in  Section  1.2  and 

,  i  A 

then  stop. 

4.  Check  Assumption  1.  If  it  holds,  determine  the  interval.  A,  where  the  Lagrange 
multiplier  resides.  If  Assumption  1  is  not  met  stop. 

5.  Apply  a  zero  finding  algorithm  to  (1.11)  within  A  to  compute  X  . 

6.  Plug  in  the  computed  X  in  (1.8)  to  determine  the  globally  optimal  solution,  x* . 
Note  that  because  I n  -  XQd  in  (1.8)  is  diagonal,  inverting  this  matrix  is 
straightforward. 

7.  Compute  the  global  minimiser  of  problem  (1.1)  using  z*  =  ST  (x*  -  Sp ) 

8.  Stop 
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6.  Localisation  of  an  Emitter  Using  TO  A  Measurements 


In  the  second  part  of  the  report  we  introduce  an  application,  which  motivated  the 
mathematical  analysis  presented  in  part  I.  The  localisation  of  an  emitter  source  or  a 
receiver  based  on  time  of  arrival  (TOA)  measurements  [8],  has  long  been  of  interest  both  in 
civilian  and  military  applications  (e.g.,  GPS).  An  objective  function  that  links  TOA 
measurements  to  the  emitter  position  in  the  plane  and  transmit  time,  t,  may  be  given  as  [9] 

i=N  2 

p(x,  y,  t)  =  Yj  ((*  -  xf  +  (y  -  T;)2  -  v2(t,  -  tf  )  (H  I) 

i= 1 


where  (x,  y)T  stands  for  the  emitter's  position  to  be  estimated  and  (x;,  y;)7  the  known  two- 
dimensional  coordinates  of  the  ith  receiver.  The  parameters  ft  are  the  measured  time  of 
arrivals.  The  parameters  <y,  are  positive  constant  weights  and  v  is  the  speed  of  light  (or 

speed  of  sound  for  an  acoustic  source).  An  estimate  of  the  emitter  position  and  transmit 
time  is  obtained  by  finding  the  global  minimum  of  the  unconstrained  cost  function  (II. 1), 
which  is  a  quartic  function  of  three  variables.  This  quartic  cost  function  is  nonconvex  and 
hence  classical  methods  such  SDP  (Semidefinite  Programming  [1][2])  based  relaxation 
algorithms  may  not  yield  the  exact  global  minimiser  and  are  generally  slower  to  execute. 

If  we  put  rt ,  —  r  =  v(tl  - 1)  then  (II.l)  can  be  written  more  compactly  as 

i=N  2 

p(x, y, r)  =  YJWi[(x- xi  f  +  ( y - yt )2  - (r  - r )2 )  (II.2) 

;=1 

If  we  expand  (II.2)  and  put 

2  e  =  x2+y2-r2  (II.3) 


in  the  polynomial  objective  function  (II.2)  then  (II.2)  reduces  to  the  four  variable  cost 
function 


i=N 


p(x,  y,  r,  e)  =  Yj  ™iP?  (*,  y,  r,  e) 

i=i 


where  pi  (x,  y,  r,e)  =  2 


r 


e-xxi-yyi  +  rri+Q- 

2  j 


(II -4) 


is  linear  in  its  four  arguments.  The 


parameter  is  defined  as  /?  =  x]  +  y2  -  r2 .  By  grouping  the  quadratic  and  linear  terms  of 

(II.2)  and  (II.4),  then  an  estimate  of  the  emitter  position  and  range  bias,  r,  can  be  obtained 
by  solving  the  equivalent  constrained  minimisation  problem 


min  zT P  z  +  2pT z  +  d 

z 

subject  to  zTQz+2qTz  =  0 


(II.5) 
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where  z  =  (x,y,r,e)T  e  R4  and  P  and  p  store  the  coefficients  of  the  quadratic  and  linear 
terms  of  (II. 4).  In  fact  as  shown  in  [9],  if  one  defines  the  ith  compound  measurement  vector 
uj  as  ui  =  2(-xj,-yi,r,  \)r  then  it  follows  that 

p  =  ^Wi(UiUf) 

i= 1 
N 

<p  =  YswiP>ui  (IL6) 

<=i 

d  =  Yj  W‘P> 

i=l 

where  d  is  a  constant  in  the  objective  function  and  therefore  can  be  ignored. 

Also  from  (II.3)  it  clearly  follows  that  Q  is  diag(l,l,-l,0)  and  q  =  (0,0,0,-l)r . 

Solving  Problem  (11.6) 

As  explained  in  Part  I  of  the  report,  problem  (II.l)  can  be  transformed  into  (1.3).  Solving  the 
optimisation  problem  (II.5)  then  amounts  to  applying  the  algorithm  of  Section  5.  Some  of 
the  steps  of  the  algorithm  call  for  checking  the  constraint  feasibility  of  (II.5)  and  whether 
the  assumptions  given  in  Part  I  are  satisfied. 

First  we  note  that  since  Q  is  indefinite  then  Qn  is  indefinite  and  hence  the  feasibility  set  is 
non-empty  and  (1.4)  is  satisfied.  Also  given  N  (  N  >  4 )  receivers  located  in  general 
positions  (e.g.,  not  accidentally  aligned),  then  P  is  positive  definite. 

Simulation  Results 

Computer  simulation  using  MATLAB  have  been  conducted  to  evaluate  the  performance  of 
the  TOA-based  geolocation.  Comparisons  were  made  with  a  similar  localisation  technique 
based  on  time-difference-of-arrival  (TDOA)  measurements.  The  TDOA-based  localisation 
method  is  selected  from  [9,  Section  3.1]  and  is  the  constrained  weighted  least  square  error 
(CWLS)  method.  This  localisation  technique  requires  solving  a  quartic  equation  in  the 
Lagrange  multiplier,  77 ,  (refer  to  (29)  in  [10]).  The  localisation  solution  is  given  by 
evaluating  (28)  in  [10].  The  TDOA  measurement  error  covariance  matrix  used  is 


“2 

1  • 

••  1 

a2 

1 

2 

1 

Q  —  TDOA 

2 

1 

1 

1  • 

••  2 

where  a  common  TOA  variance  error  is  assumed.  The  main  difference  between  the  TDOA 
and  the  TOA-based  geolocation  technique  is  the  estimation  of  the  additional  parameter, 
range  bias  (or  transmit  time). 
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The  emitter  is  at  [1,2]Km 


Figure  2:  M  ean  square  range  error  as  a  function  of  measurement  error  for  both  T OA  and  TDOA 
con  strai  n  t-  based  I  ocal  i  sati  on 

However,  it  is  emphasized  in  [11]  that  both  localisation  techniques  are  equivalent  if 
accurate  measurement  error  statistics  were  used  for  both  methods  and  that  performance 
variations  are  mainly  due  to  inaccurate  knowledge  of  measurement  error  statistics  or 
simplified  implementations.  As  can  be  seen  in  Figure  2,  the  performance  of  both 
constrained  estimation  methods  approached  the  corresponding  Cramer-Rao  lower  bound 
(CRLB)  for  a  wide  range  of  measurement  errors.  We  assume  in  this  figure  that  the  TOA 
measurement  error  is  normally  distributed  and  that  5  stationary  receivers  are  available 
and  are  positioned  respectively  at  (3,0)T  ,  (-2,1)T  ,  (0,-2)T  ,  (1,3)T  and  (3,3)T.  The  emitter  or 
source  is  located  at  (1,2)T  and  all  distances  are  expressed  in  Km.  The  errors  on  both  axes  in 
the  figure  are  displayed  in  dBm2.  In  particular  the  y-axis  of  the  figure  gives  the  mean 
square  distance  error  between  the  estimated  and  the  true  position. 


7.  Conclusions 


This  report  lays  down  a  theoretical  framework  for  minimising  a  quadratic  objective 
function  subject  to  a  quadratic  equality  constraint.  We  assume  that  the  Hessian  of  the 
objective  function  is  positive  definite.  A  detailed  algorithm  is  presented,  which  computes 
the  minimiser  in  a  straightforward  manner,  despite  the  problem  non-convexity.  Most 
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importantly  this  algorithm  finds  the  global  minimiser  exactly,  subject  only  to  computer 
hardware  precision,  and  does  not  require  special  solvers  based  on  linear  matrix 
inequalities  (LMI)  or  SDP  relaxation. 

This  work  is  motivated  by  the  need  to  solve  many  engineering  applications  such  as  time  of 
arrival  localisation.  The  localisation  of  an  emitter  source  based  on  time  of  arrival  (TOA) 
measurements,  was  the  focus  study  of  Part  II  (Section  6)  in  this  report.  It  is  shown  in  detail 
how  the  problem  is  converted  into  a  quadratic  minimisation  problem  with  one  quadratic 
equality  constraint.  The  global  solution  follows  from  the  developed  theory  of  Section  5. 
Finally  we  should  point  out  that  global  minimisation  problems  involving  more  than  one 
quadratic  constraint  are  extremely  hard  to  solve  [3]  [5]  and  until  today  no  sure  techniques 
exist  that  can  handle  more  than  one  quadratic  constraint. 
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Appendix  A: 


The  main  point  of  this  appendix  is  to  show  that  if  cmax  =  0  or  cmm  =  0  then  the  problem 
dimension  is  reduced  from  n  to  n  with  n'  <  n .  To  easily  demonstrate  this  we  assume  that 
only  cmax  =  0.  We  define  QR  to  be  the  submatrix  of  QD,  which  excludes  the  eigenvalue 


Anax  •  Clearly  the  constraint  equation  of  (1.3)  can  be  expressed  as 

xi  4  \-  Xj  —  xr  4~2cfixfl  4-5] 

1  max 


(A.1) 


where  the  index  j  (  j  >  1  )  stands  for  the  multiplicity  of  the  eigenvalue,  //max  ,  and  cR  is  the 
subvector  of  C  corresponding  to  QR.  The  vector  xR  consists  of  the  remaining  ( n  —  j ) 

elements  of  x.  Note  that  the  diagonal  matrix,  - QR ,  has  all  its  diagonal  terms  less  than  1 

Amax 

and  therefore  DR  =  I. - - —  QR  is  positive  definite.  The  objective  function,  which  can 

Amax 

now  be  expanded  as 


XT X  =  (xj2  4 - f  Xj)  4  XTrXr 

becomes  after  substitution  of  (A.l)  xRDR  xR 
minimisation  problem  reduces  to 


Ar 


s 


An 


Hence  the  new 


T  2  T  ^  s 

minimize  xrurxr  ~  crxr  ~ 

xff  A max  A max  (A.  2) 

subject  to  xtrQr  xr  +2ctrxr  +  s  <  0 


where  the  inequality  constraint  in  (A.2)  is  the  result  of  the  right  hand  side  of  (A.l)  being 
necessarily  non-negative. 


Now  if  the  unconstrained  solution  of  (A.2)  satisfies  the  constraint  xRQR  xR+2cRxR  +  s  <  0 , 


then  the  minimiser  of  (1.3)  consists  of  all  vectors 


,  where  xR  is  the  unconstrained 


solution  of  (A.2)  and  xt  •  =  (xl,---,xj)T  satisfying  (A.l).  It  clearly  follows  that  if  the  strict 

inequality  xRQR  xR+2cRxR  +  s  <  0  holds  then  multiple  solutions  exist.  If  the  unconstrained 
solution  is  not  feasible,  then  the  inequality  constraint  must  be  active  (i.e.  become  an 
equality  constraint)  and  the  optimisation  problem  then  reads  as 
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T  r\  2  T  S 

minimize  xkurxr  crxr 

X  /“max  /“max  (A.3) 

subject  to  xtrQr  xr+2 crxr  +5  =  0 

But  problem  (A.3)  is  very  similar  to  (1.1)  with  DR  positive  definite  as  required  by 
Assumption  1.  The  minimiser  of  (A.3)  is  now,  however,  sought  in  the  reduced  space  R"  7 , 
instead  of  R" . 
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